home *** CD-ROM | disk | FTP | other *** search
/ Super PC 33 / Super PC 33 (Shareware).iso / spc / sonido / timidity / source / filterin.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-08-03  |  4.4 KB  |  173 lines
  1. /*
  2.  
  3.     TiMidity -- Experimental MIDI to WAVE converter
  4.     Copyright (C) 1995 Tuukka Toivonen <titoivon@snakemail.hut.fi>
  5.  
  6.     This program is free software; you can redistribute it and/or modify
  7.     it under the terms of the GNU General Public License as published by
  8.     the Free Software Foundation; either version 2 of the License, or
  9.     (at your option) any later version.
  10.  
  11.     This program is distributed in the hope that it will be useful,
  12.     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13.     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14.     GNU General Public License for more details.
  15.  
  16.     You should have received a copy of the GNU General Public License
  17.     along with this program; if not, write to the Free Software
  18.     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  19.  
  20.    filtering.c (Author: Vincent Pagel ( pagel@loria.fr ))
  21.  
  22.    implements fir antialiasing filter : should help when setting sample
  23.    rates as low as 8Khz.
  24.  
  25.    */
  26.  
  27. #include <stdio.h>
  28. #include <string.h>
  29. #include <math.h>
  30. #include <stdlib.h>
  31.  
  32. #ifndef PI
  33.   #define PI 3.14159265358979323846
  34. #endif
  35.  
  36. #include "config.h"
  37. #include "common.h"
  38. #include "controls.h"
  39. #ifdef __WIN32__
  40. #include "instrume.h"
  41. #include "filterin.h"
  42. #else
  43. #include "instruments.h"
  44. #include "filtering.h"
  45. #endif
  46.  
  47. /*  bessel  function   */
  48. static float ino(float x)
  49. {
  50.     float y, de, e, sde;
  51.     int i;
  52.     
  53.     y = x / 2;
  54.     e = 1.0;
  55.     de = 1.0;
  56.     i = 1;
  57.     do {
  58.     de = de * y / (float) i;
  59.     sde = de * de;
  60.     e += sde;
  61.     } while (!( (e * 1.0e-08 - sde > 0) || (i++ > 25) ));
  62.     return(e);
  63. }    
  64.  
  65. /* Kaiser Window (symetric) */
  66. static void kaiser(float *w,int n,float beta)
  67. {
  68.     float xind, xi;
  69.     int i;
  70.     
  71.     xind = (2*n - 1) * (2*n - 1);
  72.     for (i =0; i<n ; i++) 
  73.     {
  74.         xi = i + 0.5;
  75.         w[i] = ino((float)(beta * sqrt((double)(1. - 4 * xi * xi / xind))))
  76.         / ino((float)beta);
  77.     }
  78. }
  79.  
  80. /*
  81.  * fir coef in g, cuttoff frequency in fc
  82.  */
  83. static void designfir(float *g , float fc)
  84. {
  85.     int i;
  86.     float xi, omega, att, beta ;
  87.     float w[ORDER2];
  88.     
  89.     for (i =0; i < ORDER2 ;i++) 
  90.     {
  91.         xi = (float) i + 0.5;
  92.         omega = PI * xi;
  93.         g[i] = sin( (double) omega * fc) / omega;
  94.     }
  95.     
  96.     att = 40.; /* attenuation  in  db */
  97.     beta = (float) exp(log((double)0.58417 * (att - 20.96)) * 0.4) + 0.07886 
  98.     * (att - 20.96);
  99.     kaiser( w, ORDER2, beta);
  100.     
  101.     /* Matrix product */
  102.     for (i =0; i < ORDER2 ; i++)
  103.     g[i] = g[i] * w[i];
  104. }
  105.  
  106. /*
  107.  * FIR filtering -> apply the filter given by coef[]
  108.  */
  109. static void filter(int16 *result,int16 *data, int32 length,float coef[])
  110. {
  111.     int32 sample,i;
  112.     int16 peak = 0;
  113.     float sum;
  114.  
  115.     for (sample = 0; sample < length - ORDER ; sample++ )
  116.     {
  117.         sum = 0.0;
  118.         for (i = 0; i < ORDER ;i++) 
  119.         sum += data[sample + i] * coef[i];
  120.         
  121.         /* Saturation ??? */
  122.         if (sum> 32767.) { sum=32767.; peak++; }
  123.         if (sum< -32768.) { sum=-32768; peak++; }
  124.         result[sample] = (int16) sum;
  125.     }
  126.     
  127.     if (peak)
  128.     ctl->cmsg(CMSG_ERROR, VERB_NORMAL, 
  129.           "Saturation %2.3f %%.", 100.0*peak/ (float) length);
  130. }
  131.  
  132. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  133. /* Prevent aliasing by filtering any freq above the output_rate        */
  134. /*                                                                     */
  135. /* I don't worry about looping point -> they will remain soft if they  */
  136. /* were already                                                        */
  137. /*_____________________________________________________________________*/
  138. void antialiasing(Sample *sp, int32 output_rate )
  139. {
  140.     int16 *temp;
  141.     int i;
  142.     float fir_symetric[ORDER];
  143.     float fir_coef[ORDER2];
  144.     float freq_cut;  /* cutoff frequency [0..1.0] FREQ_CUT/SAMP_FREQ*/
  145.  
  146.  
  147.     ctl->cmsg(CMSG_INFO, VERB_NOISY, "Antialiasing: Fsample=%iKHz",
  148.           sp->sample_rate);
  149.  
  150.     /* No oversampling  */
  151.     if (output_rate>=sp->sample_rate)
  152.     return;
  153.     
  154.     freq_cut= (float) output_rate / (float) sp->sample_rate;
  155.     ctl->cmsg(CMSG_INFO, VERB_NOISY, "Antialiasing: cutoff=%f%%",
  156.           freq_cut*100.);
  157.  
  158.     designfir(fir_coef,freq_cut);
  159.     
  160.     /* Make the filter symetric */
  161.     for (i = 0 ; i<ORDER2 ;i++) 
  162.     fir_symetric[ORDER-1 - i] = fir_symetric[i] = fir_coef[ORDER2-1 - i];
  163.     
  164.     /* We apply the filter we have designed on a copy of the patch */
  165.     temp = safe_malloc(sp->data_length);
  166.     memcpy(temp,sp->data,sp->data_length);
  167.     
  168.     filter((int16 *)sp->data,temp,sp->data_length/sizeof(int16),fir_symetric);
  169.     
  170.     free(temp);
  171. }        
  172.  
  173.